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Introduction 


In  recent  years  a  lot  of  attention  has  been  paid  to  the  analysis  of  the  structures  and  parts 
operating  under  extreme  thermal  and  mechanical  loads.  Special  interest  has  arisen  in  the 
high  temperature  turbine  parts  which  subject  to  long  term  viscoplastic  deformations  as 
well  as  time  independent  inelasticity  caused  by  high  level  of  mechanical  stress.  This 
combination  of  creep  and  plasticity  leads  to  the  damage  nucleation  and  growth  and  a 
significant  reduction  of  the  part  expected  service  life.  Accurate  prediction  of  material 
response  for  combined  creep  and  plasticity  defonnations  is  a  complex  problem  greatly 
compounded  during  cyclic  loading.  During  the  air  foil  service  the  strain  controlled  cyclic 
elastic-plastic  deformations  corresponding  to  airplane  maneuvers  are  combined  with  so- 
called  dwell  times  causing  viscoplastic  effects  such  as  creep  and/or  stress  relaxation.  The 
development  of  unified  creep  -  plasticity  model  capable  of  predicting  cyclic  non- 
isothermal  loading  conditions  is  of  extreme  importance. 

Our  constitutive  models  are  based  on  state  variables  representing  creep,  plasticity  and 
damage  evolution.  The  basic  physical  mechanism  of  the  inelastic  deformation  analyzed  in 
this  paper  is  dislocation  glide  or  so-called  crystallographic  slip.  We  study  the  Ni-based 
superalloy  PWA  1484  having  LI 2  crystallographic  structure.  Deformation  occurs  along 
12  octahedral  slip  systems  {111}<110>  and  along  6  cube  slip  systems  {00 1  }<1 10>. 

Since  the  cube  systems  slip  resistance  is  higher  than  the  octahedral  systems,  the  cube 
systems  activity  is  limited.  However,  at  high  temperature  and  for  initial  crystallographic 
orientations  close  to  <1 1 1>  the  y’  particles  are  sheared  by  cube  systems.  (C.  Allan  1995, 
Stouffer  and  Dame  1996).  Please  see  Table  1  for  the  slip  system  nomenclature  in  a  single 
crystal  superalloy. 

Because  of  the  wide  temperature  range  the  parts  have  been  exposed  during  the  service  the 
inelastic  defonnation  might  be  decomposed  into  creep  and  plasticity  parts.  At  high 
temperatures  the  creep  is  a  dominant  mode  of  deformation  and  due  to  stress  re¬ 
distribution  the  applied  stress  level  remains  at  moderate  levels  and  does  not  cause  plastic 
deformation.  The  description  of  the  creep  model  is  given  in  Staroselsky  and  Cassenti 
(2008).  However,  when  the  part  is  being  cooled  down  creep  simply  disappears  and 
inelastic  deformation  takes  place  by  crystallographic  slip  which  does  not  vary  with 
temperature  as  strongly  as  creep.  Our  creep  experiments  show  that  creep  rates  follow 
Arrhenius  type  relationship  with  the  value  of  the  activation  energy  varying  in  the  range 

from  (9<001>  =  6.97 E  - 19  to  Q<U1>  =  7.30 E  - 19  ^/atom  with  some  stabilization 

(decrease  of  the  activation  energy)  for  the  temperatures  below  700  C.  Plasticity  plays  the 
major  role  in  residual  stress  generation  during  the  cooling  stage  of  TMF  process. 

In  a  rate  dependent  approach  a  constitutive  law  for  the  slip  rate  is  assumed,  usually  in  the 
power  law  form  (for  example,  Asaro  R.J.  1983,  Asaro  R.J.  ,  Needleman,  A.  1985, 
Kalidindi  et  al.  1992).  This  means  that  slip  systems  are  active  even  before  the  resolved 
shear  stress  reaches  critical  values.  The  rate  of  inelastic  defonnation  is  proportional  to  the 
power  of  the  effective  resolved  shear  stress  and  the  corresponding  boundary  value 
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problem  has  a  unique  solution.  The  higher  the  exponent  in  the  power  law  the  closer  this 
approach  gets  to  the  classical  plasticity  theory.  However,  the  price  one  has  to  pay  for  this 
is  the  increase  in  computational  stiffness  of  the  problem  and  an  increase  in  the 
computational  time.  Nevertheless  rate-dependent  power  law  has  demonstrated  its 
efficiency  and  has  been  widely  applied  for  both  monotonic  and  cyclic  problems  of 
plasticity.  However,  when  plastic  analysis  is  combined  with  creep  predictions,  the  rate 
dependent  nature  of  the  plastic  part  of  the  flow  rule  might  cause  some  problems.  A  very 
good  example  is  the  inelastic  analysis  of  single  crystals  of  Ni-based  superalloys  which 
has  a  high  creep  rate  at  stresses  of  the  order  of  two  thirds  of  the  yield.  Let  us  consider  a 
typical  (Nissley,  D.  et  al  1991,  Staroselsky  and  Cassenti,  2008)  rate  depended  model 
where  the  set  of  calibration  parameters  are  as  follows: 


r  =r  o 


,  where  ;/(l  =  0.0001  sec  1  and  m=30.  It  is  important  to  note  that  there  are 


y 

only  two  independent  parameters  in  such  a  formulation:  —  and  the  exponent  m. 

sm 

Calibration  is  based  on  the  computational  predictions  of  plastic  response  to  be  in  a  good 
accord  with  experimental  observations.  Then,  after  100  hours  of  creep  at  the  stress  level 
of  80%  of  yield  the  rate  depended  plastic  tenn  will  predict  the  “creep”  of  approximately 
5%  ,  which  is  more  than  a  quarter  of  total  observed  inelastic  deformation  to  failure.  The 
plastic  term  predictions  interfere  with  creep  prediction  the  closer  the  applied  load  is  to 
yield.  This  problem  might  be  eliminated  by  increasing  the  exponent  to  the  level  of  80 
and  higher,  as  shown  in  (Kalidindi  and  Anand,  1994)  but  this  leads  to  its  greater 
computational  difficulties.  Therefore  there  is  a  need  for  the  further  development  of  rate- 
independent  crystal  plasticity  and  combination  with  viscous  terms.  Also,  it  is  important 
to  note  that  rate  sensitivity  of  Ni-based  superalloys  is  small,  which  in  turn  leads  to  high 
values  of  the  exponent. 


We  use  two  approaches  to  describe  the  rate  independent  behavior  of  single  crystal 
materials.  The  first  one  is  based  on  the  concept  of  the  potentially  active  systems  and  is 
reduced  to  the  system  of  linear  equations  for  the  slip  shear  increments.  The  trial  elastic 
increment  is  used  to  select  the  potentially  active  slip  systems  and  to  construct  the  linear 
equations.  The  idea  of  the  method  was  originally  formulated  for  the  slip  systems  in  FCC 
crystals  under  monotonic  loading  conditions  by  Anand  and  Kothari  (1996)  and 
generalized  for  combination  of  slip  and  twin  systems  by  Staroselsky  and  Anand  (1998). 
The  second,  more  general  approach  uses  the  idea  that  the  plastic  strain  rate  is 
proportional  to  the  total  energy  rate.  Basically,  by  cancelling  the  derivatives,  we  come  to 
the  rate-independent  formulation.  Some  of  the  earliest  possible  formulations  were 
proposed  by  Kremple  (1976),  Valanis  (1976),  and  Walker  (1980).  Later  an  isotropic 
version  developed  by  Cassenti  (1983)  was  shown  to  reduce  in  the  limit  to  classical 
plasticity.  The  formulation  required  only  two  new  temperature  dependent  material 
parameters.  One  was  the  yield  stress  and  the  other  was  a  curvature  parameter  that  can 
provide  a  gradual  transition  form  a  smooth  uniaxial  stress-strain  curve  to  a  stress-strain 
curve  with  a  sudden  change  from  the  elastic  slope  to  a  perfectly  plastic  response.  Both 
isotropic  and  kinematic  hardening  were  included  using  the  time  dependent  creep 
formulation.  We  generalize  the  formulation  of  Cassenti  to  single  crystal  materials  and 
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add  an  appropriate  power  law  to  make  the  formulation  agree  exactly  with  a  perfectly 
plastic  formulation  described  in  the  first  method.  In  this  paper  we  focus  on  the 
development  of  generic  rate-independent  crystal  plasticity  models  suitable  for  cyclic  non- 
isothermal  elastic-plastic-viscoplastic  analysis. 

The  plan  of  the  paper  is  as  follows:  We  start  section  2  with  generalization  of  the  rate- 
independent  crystal  plasticity  model  developed  by  Anand  and  Kothari  (1996)  to  cyclic 
deformations.  Then,  we  formulate  the  new  rate-independent  model  which  combines  the 
advantages  of  both  rate-dependent  power  method  and  with  a  rate-independent  approach. 
We  will  compare  the  computational  results  obtained  by  all  of  these  approaches  for  corner 
(symmetric)  and  a  non-corner  (single  slip)  initial  crystal  orientations.  We  will  discuss  slip 
system  activities  for  different  crystallographic  orientations  including  octahedral  and  cube 
slip  systems  in  the  section  4.  In  the  next  section  we  will  apply  the  developed  rate- 
independent  plasticity  approach  to  the  modeling  of  simple  tension  specimen  and  compare 
texture  evolution,  numerical  stability,  and  necking  effects  with  homogeneous  numerical 
solutions  and  available  experimental  results.  Section  6  provides  results  of  model 
predictions  of  cyclic  deformations  of  Ni-based  superalloy  and  compares  them  against 
coupon  testing  results.  We  close  the  paper  with  some  concluding  remarks. 


2.  Governing  Equations 

The  overall  plastic  response  is  taken  as  a  sum  of  responses  from  small  regions  of  a  single 
crystal  playing  the  role  of  representative  volume  elements  (RVE).  The  deformation  of  a 
crystal  is  taken  as  the  sum  of  contributions  from  overall  elastic  distortion  and  plastic 
deformation. 

The  governing  variables  in  the  constitutive  model  are  as  follows:  the  Cauchy  stress 
tensor  T  and  the  defonnation  gradient  F  .  Each  crystal  slip  system  is  specified  by  a  unit 
normal  to  the  slip  plane  and  a  unit  vector  along  the  slip  direction  . 

The  total  deformation  gradient  is  multiplicatively  decomposed  on  elastic  and  plastic  parts 
as 

F=FeFp  and  detFp=l.  (1) 

The  elastic  deformation  gradient  F%  (detFe  >  0)  describes  elastic  distortion  and  gives 
rise  to  the  stress  T  .  The  constitutive  equation  for  the  second  Piola-Kirchoff  stress  tensor 
is  taken  as  a  linear  relation: 

T*  =  c[Ee]  where  Ec(r)  =  i(ce(r)  -i);  and  Ce(r)  =  FeTFe  (2) 

at  the  current  moment  of  time  denoted,  x. 

The  Cauchy  stress  tensor  is  the  work  conjugate  stress  corresponding  to  the  Cauchy  - 
Green  elastic  strain  measure  and  is  calculated  as  follows: 


1  Generally  speaking  the  elastic  relation  looks  as  follows:  T* 


L  Ee  —  A(0  —  0O  )j  where  A  is  thermal 


expansion  tensor,  ©  is  current  temperature  and  0O  is  the  reference  temperature.  All  results  we  report  in  this 
paper  are  valid  and  derived  general  formulae  and  algorithms  do  not  change  when  we  add  this  thermal  strain 
term.  For  the  sake  of  simplicity  we  will  use  formula  (2)  in  this  paper. 
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T  = 


1 


-FeT'Fe 


_  (3) 

det(Fe) 

The  evolution  equation  for  the  viscoplastic  defonnation  gradient  is  given  by  the  flow 
rule: 


ip  =  I/Fp  where  If  = 

slip  systems 


and  S“  =  m“ 


n 


(4) 


The  shear  rate  along  each  slip  system  ya  is  given  in  terms  of  the  resolved  shear  stress 

(RSS)  r  =  T*  •  S“  ,  slip  systems  resistances  and  equilibrium  back  stress.  Evolution  of 
crystallographic  texture  is  explicitly  defined  by  the  elastic  part  of  defonnation  gradient. 


a 

m,  =  r  m„ 

a  T^e-1  a 

nt  =F  n0 


(5) 


In  order  to  complete  the  calculations  one  has  to  be  able  to  calculate  slip  systems  shear 
rate  or  slip  systems  shear  increment  in  an  incremental  formulation.  In  rate-dependent 
model  this  problem  has  been  resolved  by  defining  a  power  law  viscoplastic  relation  with 
a  back  stress:  for  the  inelastic  strain  rate  along  a-th  slip  system: 


=  fo 


(1  ~dv) 


-co 


sgn  (ra-coa). 


(6) 


where  y0  is  a  temperature  dependent  time  constant,  and  sa  is  the  deformation  resistance 

OC 

of  a-th  slip  system;  p  is  resolved  shear  stress,  of  is  the  slip  system  back  stress,  n  is 
the  creep  exponent  which  is  very  high  (~80)  to  simulate  plastic  response,  dv  is  the  void 
related  damage  parameter,  and( ' )  is  the  rate  of  change  with  respect  to  time. 

Next,  latent  hardening  evolution  has  been  described  by  modifying  Asaro’s  (1983) 

f  a\P 

hardening  rule  sa  =  h0 


i-^ 


rf  ,  with  hardening  matrix 


h"/!  =  \q  +  (1  -  q)8af> }  for  temperature  dependent  h0  and  s  * .  We  model  cyclic  effects  by 
defining  for  each  slip  system  a  specific  internal  equilibrium,  or  back  stress  co.  The  back 

Q 

stress  has  a  limiting  saturation  value  com  =  — ,  corresponding  to  the  stabilization  of  the 

C2 

equilibrium  portion  of  the  stress  (i.e.,  the  portion  of  the  stress  that  has  not  yet  contributed 
to  the  inelastic  strain)  and  evolves  according  to  the  following  relationship  (Nissley,  D.,  et 
al,  1991,  Voyiadjis,  G.Z.,  Huang,  W.,  (1996),  Stouffer,  D.C.,  and  Dame,  L.T  (1996)): 

coa  =cja -c2\ya\coa  (7) 

Equation  (7)  requires  two  additional  coefficients  c,  and  c,  that  are  explicit  functions  of 
temperature.  The  effect  of  voids  on  the  slip  rate  can  be  calculated  through  a  correction  of 
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the  effective  stress: 


(!-<*„) 


-  ( o 


(an  effective  increase  in  the  resolved  shear  stress).  This 


is  equivalent  to  a  Kachanov-type  damage  parameter  (Kachanov  1986,  Lemaitre  1996)  . 
For  the  sake  of  simplicity  we  will  not  consider  the  effects  of  voids  in  this  paper. 


3.  Generalization  of  incremental  linear  system  formulation  to  cyclic 
plasticity. 

One  of  the  major  problems  of  crystal  plasticity  is  to  determine  the  unique  set  of  the  active 
slip  systems  and  their  corresponding  shear  rates  at  each  time  increment.  At  each  strain 
increment,  only  five  strain  components  can  be  uniquely  prescribed.  Because  the  number 
of  slip  systems  for  FCC  systems  are  12  and  for  LI  2  structures  are  18,  the  problem  of 
selecting  active  ones  is  crucial.  In  a  power  law  model  all  slip  systems  are  active  just  the 
value  of  slip  rate  might  be  infinitesimal  if  the  resolved  shear  stress  is  small.  If  the  value 
of  exponent  is  small  (from  3  to  5),  which  is  the  case  for  creep  phenomena,  then  the  slip 
systems  with  significant  but  not  the  largest  values  of  the  resolved  shear  stress  (  RSS)  are 
still  active  and  the  total  shear  along  these  “second  order”  slip  systems  could  be 
significant.  In  the  rate  independent  (RI)  model,  the  system  might  be  active  only  if  the 
effective  RSS  is  equal  to  slip  resistance  (i.e.,  the  yield  stress).  Then 

ya  >  0  if  and  only  if  |r"  -  coa  |  =  sa  (r) .  (8) 

In  this  section  we  will  reduce  the  problem  to  the  systems  of  the  linear  equations  with 
respect  to  shear  increments.  This  system  should  be  solved  at  each  incremental  step  using 
any  numerical  method  minimizing  errors  for  example,  the  method  of  singular  value 
decomposition  (SVD),  which  will  define  the  unique  solution  for  the  plastic  strain  in  L2 
norm  sense  (Anand  and  Kothari  1996)  or  penalty  functions  (Staroselsky  and  Anand 
1998). 

We  use  an  incremental  formulation  where  we  define  all  parameters  at  the  end  of  the  time 
increment  r  =  t  +  At  based  on  state  in  the  beginning  of  the  increment  at  time  t. 

Expression  (4)  has  the  incremental  form  as  follows: 

Fp(t)  =  V’  (r)F  '  (0  =  (1  +  X  Ar ' X  )FP  (0  (9) 

a 

We  need  to  express  elastic  Cauchy-Green  tensor  through  total  and  plastic  deformation 
gradients  using  equations  (1)  and  (2)  as: 

Ce(r)  =  F^TFrFF^1  (10) 

Numerical  integration  is  taken  using  a  forward  Euler  integration  scheme.  In  order  to 
perform  it  within  this  RI  model,  we  define  a  trial  step  where  we  assume  that  the  whole 
increment  is  elastic.  The  stress  overshoot  is  compensated  for  by  crystallographic  slip.  The 
constitutive  relations  for  the  trial  step  are  taken  as  follows: 
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(11) 


T*"'  =  z[Ee"  (r)], 

Ee  tr(r)  =  y(ce,r(r)  -i),  and 


Where 

Fe"(r)  =  F(r)Fp(ZL)  1  (12) 

and  the  trial  resolved  shear  stress  along  each  slip  system  is 

Tatr  =  T  *tr  *S“  .  (13) 

Substituting  (12)  into  (11)  and  using  expressions  (9)  and  (10)  one  can  obtain  the 
relationship  between  true  and  trial  stress  tensors  (Kalidindi,  et  al  1992),  which  already 
accounts  for  the  amount  of  incremental  shear,  as: 

x*  =  x **  L,\sym(cetr Sq  %ya  (14) 

a 

and  similarly  for  the  resolved  shear  stresses 

=  Tatr  ~Y^ym(ce,rslh:Ayfj  = 

=  Tair  -  ^Z,(jy/«(cefrSS  ))S“ \Ay\p signer p,r  -a>fi)  ^ 

P 

At  the  next  step,  following  the  idea  of  Anand  &  Kothari  (1996)  we  assume  that  during 
small  strain  increment,  changes  in  the  slip  directions  defined  as 
sign(Ay“  (r))  =  sign(r“  (r)  -  of  (r))  would: 

(i)  be  the  same  for  the  true  and  trial  increments 
sign(r“  -  of  )  =  sign(r“"'  -  oftr  \  and 

(ii)  the  effect  of  change  in  back  stress  values  during  one  increment  does  not  affect 
the  slip  direction:  sign(r“(r)  -  coa {z))=  sign(r“(t)  -  coa (t)). 

We  need  these  assumptions  in  order  to  compute  absolute  values  of  the  effective  resolved 
shear  stress  used  for  the  yield  criteria 

\ra -a>a\  =  sa(r)  (8) 

The  systems  for  which  trial  values  of  effective  shear  stress  are  less  than  slip  resistance  in 
the  beginning  of  increment  s(t)  are  called  inactive  and  are  not  analyzed  during  the 
increment.  For  the  other  or  so-called  “potentially  active”  slip  systems,  the  absolute  value 
of  the  effective  shear  stress  can  be  calculated  through  the  resolved  shear  stress  trial  value. 

|r“ -o)a\  =  \ratr  - coa (r)| - ^z(jym(ce<rSjJ ))sjj \Ay\^ sign{rp  - (Qp) sign{ra  -coa)  (16) 

p 

By  the  use  of  the  incremental  formulation  for  the  back  stress  (7): 
of  (r)  =  of  (0  +  (c,sign(Aya)-  C2coa  {t)\Aya  | 
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We  immediately  obtain  that 


\za,r  -  coa  (r)\  =  \ra,r  -coa(t)\-(ci  -C2(oa  {t)\ign(ratr  -  a*  (tj\Aya\ 
And  finally,  using  the  latent  hardening  evolution  equation  sa  (r)  =  sa 


(17) 


{t)  +  Y,hap  Ay1 


p 

we  reduce  the  problem  of  finding  absolute  values  of  slip  increments  determination  to  a 
solution  of  the  linear  system, 

B|Ay|  =  b  (18) 


where 

Bafi  =hap  +(ct  -C2(oa{t)\ign(Ta,r  -  cop) sign(zatr  -  ma) 

ba  =  \ratr -o)a{t\-sa{t).  (19) 


This  system  has  a  number  of  variables  equal  to  the  number  of  “potentially  active”  slip 
systems,  however  not  all  of  the  linear  equations  are  always  linearly  independent.  The 
system  might  be  solved  using  any  appropriate  approximation  scheme  from  linear  algebra. 
One  of  the  “proven”  and  effective  algorithms  is  SVD  which  gives  the  least  square 
solution.  We  do  not  discuss  the  advantages  and  disadvantages  of  different  linear  algebra 
algorithms  and  use  that  numerical  approach  throughout  this  paper. 


4.  Time-independent  formulation 


In  the  previous  section  we  have  described  the  rate  independent  plasticity  calculations  by 
using  approximate  solution  of  a  linear  system.  This  is  the  limiting  case  when  the  yield 
surface  is  well-defined.  However,  there  is  still  a  need  to  combine  the  advantages  of  a 
continuous  rate  independent  plastic  behavior  description.  This  is  especially  true  for 
cyclic  non-isothermal  conditions  present  in  thennal  mechanical  fatigue  (TMF).  Such  a 
formulation  would  automatically  include  the  characteristic  variation  in  yield  surfaces 
between  materials.  In  this  section  we  fonnulate  our  new  rate-independent  flow  rule  and 
compare  the  deformation  and  slip  activities  model  predictions  against  classical  power  law 
as  well  as  the  linear  system  rate-independent  method  described  above. 


We  include  a  time  independent  inelastic  strain  by  an  additive  decomposition  of  the 
deformation  velocity  gradient  into  three  parts:  elastic,  the  time  dependent  creep  strain 
rate,  and  a  time  independent  plastic  one.  We  can  insure  that  the  plastic  strain  rate  is  time 
independent  by  making  it  proportional  to  another  rate.  Let  us  consider  the  example,  a 
plastic  strain  rate  that  is  proportional  to  the  plastic  work  rate,  Wp  ,  then 


p 

>j 


dt 


f{wp)vp  =f{wp) 


dW 

dt 


p 


(20) 


and  the  differentials  will  cancel  each  other  in  an  integration  to  yield 
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(21) 


del 

dWp 


Equation  (21)  makes  £p  a  function  of  Wp  independent  of  any  rates,  and  subsequently 

the  model  is  rate-independent.  Now,  let  us  apply  this  concept  to  fonnulate  a  crystal 
plasticity  rate-independent  flow  rule. 


For  each  of  the  slip  systems  the  inelastic  strain  rate,  y  “ ,  is  a  function  of  the  effective 
stress  for  that  slip  system,  Ta-coa,  and  we  could  choose, 


ya  = - - — <WP  >  (22) 

where  (xj  is  the  unit  ramp  function  of  x  and  crr  is  a  temperature  dependent  material 

parameter.  Expression  (22)  is  deficient  since  a  slip  system  will  be  active  with  noticeable 
shear  rate  even  for  relatively  moderate  values  of  the  effective  resolved  shear  stress,  which 
is  much  less  than  one  for  actually  active  slip  systems.  This  can  be  corrected  by 
modifying  equation  (22)  through  the  use  of  a  power  law  in  the  effective  resolved  shear 
stress  as  follows: 


r 


(23) 


Here  we  used  the  fact  that  power  function  with  the  argument  varying  from  zero  to  unity 
behaves  almost  like  a  switch  function  for  high  powers  where  the  closer  to  ideal  switch 
(Heaviside  function)  the  higher  the  exponent  is.  In  our  case  we  use  values  of  the  exponent 
m  <  30  which  guarantees  good  selection  of  slip  activities  and  free  of  the  computational 
problems  specific  for  power  models  with  very  high  values  of  the  exponent  values. 


Equation  (23)  has  assumed  that  only  the  plastic  work  contributes  to  the  plastic  strain,  but 
the  total  work  rate  may  also  be  an  important  variable.  Since  Wp  =W  -  We ,  we  can  use  a 
weighted  sum  of  the  total  and  elastic  work  rates  by  replacing  the  plastic  work  rate  with  a 
weighted  sum  in  equation  (23)  in  the  following  way: 


r  = 


_  a 

a 

T  - CO 

T  -CO 

2 

^00 

S 

W-k{w~Wp )' 


(24) 


where  k  is  a  temperature  dependent  material  parameter.  The  total  work  rate  is  given  by 

^  =  Tjailv  =C7lii’  (25) 


and  the  plastic  work  rate  is 

»"  =  IE 


(26) 


Substituting  equation  (26)  into  equation  (24)  for  the  plastic  work  rate  and 
solving  equation  for  Wp  yields 
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(27) 


MX*>I 


IF 


p  _ 


r  — 


v  y 


r  —  <y 


i-*Z 


f  a  A 

2 

a 

m 

T  -CO 

T  -CO 

_ 2 

v  ) 

s 

Substituting  equation  (28)  into  equation  (24)  yields 


Y  = 


(l  k)(w\ T  ^ 

a 

T  —CO 

m 

l1  K)yy  /  9 

s 

i-*E 

a 

(  f  ~(Oa  V 

_  a  a 

T  —CO 

m 

l  j 

s 

(28) 


Identifying  cr^  with  s  in  the  equation  (28)  and  assuming  that  state  variable  s  might  be 
different  for  different  slip  systems  we  get  the  final  form  for  the  rate  independent  plastic 
strain  rate,  we  can  write 

'w' 


MX 


r 


-a  a 

T  —CO 


m+ 1 


1-*Z 


a 

T  —CO 


m+2 


(29) 


Thus,  we  formulated  the  rate-independent  flow  rule  and  introduced  one  more  parameter 
0  <  k  <  1  .  If  the  value  of  this  parameter  is  close  to  zero,  yield  values  change  gradually 
as  it  usually  observed  during  the  defonnation  hardening.  In  the  limiting  case  of  k  — >  1  (in 
Eq.  29)  there  will  be  a  sudden  change  in  the  plastic  strain  rate  whenr"  -  of  =  sa ,  and, 
hence,  5“  can  be  interpreted  as  a  yield  stress.  Interestingly  enough  that  graduate 
hardening  behavior  is  correlated  to  the  partial  amount  of  dissipated  energy  proportional  to 
the  yield  strain  increment.  Typical  stress-strain  curves  with  different  values  of  parameter 
k  are  shown  in  Fig.  1. 


5.  Model  Implementation  and  Evaluation  of  the  Parameters 

In  this  section  we  compare  the  predictions  obtained  by  our  new  rate  independent 
plasticity  model  against  predictions  obtained  by  a  power  law,  and  the  solution  of  the 
linear  system  described  in  section  2,  with  available  experimental  data.  We  discuss  single 
crystal  Ni-based  superalloy  constitutive  behavior  at  relatively  high  temperatures  where 
the  loading  applied  along  one  of  the  comer  orientation  <00 1>  or  <1 1 1>  and  along  an 
internal  orientation  of  <  236  >  and  <  1  23  > . 

Anisotropic  elastic  properties  are  defined  by  three  different  moduli  that  are  changed  with 
temperature  and  shown  in  Fig.  2.  The  fourth  order  elasticity  tensor  has  the  form  specific 
for  orthotropic  materials  with  cubic  symmetry: 
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Cm=c  11;  (1=1,3) 
C,W=C\2-  (z,  j  =  1,  3) 
C,,  =  C44;  O',  j  =  1,3) 


We  consider  single  crystal  Asaro  hardening  behavior  sa  =hG 


f 

1 - : 


V 


y 


where 

P 


the  parameters  h0  and  initial  value  of  slip  resistance  s0  for  both  octahedral  and  cube 


systems  are  also  functions  of  temperature.  Results  of  slip  resistance  calibration  against 
experimental  data  are  shown  in  Fig  3.  Due  to  low  strain  hardening  of  PWA1484  for  some 
example  calculations  we  use  a  non-hardening  model. 


Using  the  developed  constitutive  models  we  perform  simulations  of  displacement 
controlled  simple  tension  tests.  We  use  ANSYS  ET  185  8-node  brick  elements  for  both 
single  element  and  3744-element  cylinder  simulations.  Single  element  simulations  satisfy 
homogeneous  uniaxial  stress  and  strain  distribution.  Analysis  shows  that  two  types  of 
boundary  conditions  provide  the  required  uniaxial  stress-strain  distribution.  We  will  call 
the  single  element  with  the  following  boundary  conditions:  upper  and  lower  surfaces 
remain  parallel  and  orthogonal  to  the  displacement  direction,  the  element  type  A.  The 
element  with  two  vertical  sides  parallel  to  the  loading  direction  will  be  called  the  element 
type  B.  Both  elements  and  sketches  illustrating  the  deformation  of  each  type  are  shown  in 
Fig.  4  (from  Dieter,  1988). 


There  is  no  numerical  difference  between  predictions  of  single  elements  for  A  or  B  types 
subject  to  simple  tension  of  high  symmetry  corner  orientation  <00 1>  and  <1 1 1>  due  to 
highly  homogeneous  deformation.  Moreover,  the  single  element  prediction  is  exactly  the 
same  as  the  numerical  predictions  obtained  from  cylinder  calculations.  Crystal  orientation 
does  not  change  during  the  numerical  tests  in  both  corner  orientations.  There  are  eight 
equally  active  octahedral  systems  for  <00 1>  crystallographic  orientation  tests.  The 
normalized  stress-strain  curves  obtained  by  methods  described  above  and  also 
experimentally  measured  for  980  C  are  shown  in  Fig  5a.  There  are  six  equally  active 
octahedral  slip  systems  (T 1 1)[T01] ,  (1  ll)[0lT],  (TlT)[01 1] ,  (1 1 1)[10l] ,  (I  I l)[0l  I] 
and  (1  1  1)[1  10]  as  well  as  3  identically  active  cube  slip  systems  (100)[01 1] , 

(0 10)[1 0  1],  and  (00 1)[1  1  0]  observed  during  simulations  of  simple  tension  along  a 
<1 1 1>  crystallographic  orientation.  State  variables  have  been  calibrated  in  such  a  way 
that  cube  slip  system  shear  is  two  orders  of  magnitude  larger  than  octahedral  slip  shear. 
The  normalized  stress-strain  curves  for  this  crystallographic  orientation  at  980  C  are 
shown  in  Fig  5b. 


When  loaded  along  a  central  direction  inside  the  primary  stereographic  triangle  the  single 
crystal  deforms  initially  along  a  single  slip  system  with  the  highest  value  of  resolved 
shear  stress  (i.e.,  with  the  highest  value  of  Schmidt  factor).  Usually  such  a  deformation  is 
very  unstable.  We  conducted  simple  tension  simulations  (type  A)  of  a  single  crystal  along 
<  236  >  crystallographic  direction.  The  stress-strain  curves  for 
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980  C  are  shown  in  Fig  6.  The  slip  system  activity  is  shown  in  Fig.  7.  Deformation  starts 
with  a  single  octahedral  slip  system (1 1 1)[  1  01],  then  when  the  strain  reaches  about  8% 
the  second  conjugate  slip  system  (111  )[01 1]  becomes  active.  It  is  important  to  note  that 
the  cube  slip  system  (00 1)[1  1  0]  is  also  active  from  the  beginning,  as  can  be  seen  from 
the  plots  in  Fig.  7,  but  the  total  amount  of  cube  slip  system  shear  is  much  less  than  due  to 
octahedral  slip  systems.  The  slip  resistance  of  the  conjugate  slip  system  is  bigger  than  the 
slip  resistance  of  the  active  one  due  to  latent  hardening.  The  effective  resolved  shear 
stress  reaches  the  critical  value  at  deformation  levels  slightly  larger  than  for  the  non¬ 
hardening  case'  and  subsequently  the  re-orientation  slightly  overshoot  the  stereographic 
triangle  boundary  [001]  -  [  1  11].  After  activation  of  the  conjugate  slip  system  (1 1  1  )[0 1 1] 
the  lattice  rotates  toward  the  [111]  pole  along  the  triangle  boundary.  The  prediction  of 
our  rate-independent  model  is  indistinguishably  close  to  the  predictions  obtained  by  the 
solution  of  the  linear  system  and  by  the  power  law  (rate-dependent  model)  with  the 
exponent  higher  than  80.  It  is  important  to  note  that  the  rate-dependent  power  law  with 
small  values  of  exponent  specific  for  creep  (3  to  5)  predicts  that  more  slip  systems  are 
active  leading  to  more  homogeneous  deformation  and  slower  crystal  re-orientation.  For 
example,  for  creep  along  <  236  >  crystallographic  orientation  slip  systems  (1 1 1)[01  1] 
and  (11  1 )[  1  0  1  ]  are  also  active,  however  total  slip  along  these  systems  is  significantly 
smaller  than  along  the  primary  and  the  conjugate  ones. 

Deformation  along  <  236  >  is  very  sensitive  to  the  boundary  conditions.  The  stress- 
strain  relationship  obtained  under  different  types  of  boundary  conditions,  such  as  Type  A 
as  was  reported  above,  type  B  loading  as  well  as  numerical  results  from  the  cylinder 
specimen  model  is  shown  in  Fig.  8.  Deformation  takes  place  along  the  single 
(1 1 1)[  101]  slip  system  till  collapse  occurs  by  excessive  deformation.  Crystal 
reorientation  is  rotated  slightly  toward  the  [001]  pole  from  the  boundary  conditions  as 
shown  in  Fig.  9  together  with  prediction  obtained  for  type  A  loading.  We  evaluated  the 
one  element  prediction  by  comparing  the  numerical  results  against  our  3744-element 
cylinder  specimen  model.  The  top  and  bottom  of  the  numerical  specimen  are  rigidly 
prevented  from  rotation  and  absolutely  symmetric.  In  order  to  eliminate  end-effects  the 
last  2  layers  of  finite  elements  from  each  end  are  made  elastic.  We  use  35  elastic-plastic 
central  layers  of  finite  elements  and  96  elements  in  each  cross-section  as  shown  in  Fig  10. 
Deformation  is  very  sensitive  to  a  minor  variation  of  boundary  conditions  and  starts 
necking  at  deformation  levels  of  about  7-10%  as  shown  in  Fig  11.  An  increase  in  the 
number  of  elements  at  least  within  the  same  order  of  magnitudes  does  not  change  the 
character  of  deformation  which  appears  to  be  very  inhomogeneous.  The  central  elements 
deform  almost  according  to  loading  scheme  A,  while  surface  and  near  surface  elements 
clearly  deforms  as  type  B  elements.  Due  to  inhomogeneities  in  the  deformation  the 
crystal  lattice  reorientation  of  the  central  element  is  different  than  type  A  results  but 
reasonably  close  to  it  as  can  be  seen  from  Fig.  9.  Due  to  the  lack  of  the  symmetry  during 
the  deformation  the  model  predicts  noticeable  ovalization  of  the  cross-section  as  shown 


'  We  performed  modeling  for  the  “theoretical”  non-hardening  case  and  the  conjugate  slip  system  became 
active  exactly  on  the  boundary  of  stereographic  triangle. 
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in  Fig.  12.  Such  an  ovalization  is  specific  for  any  non-comer  crystallographic  orientation 
and  might  be  used  for  specimen  life  prediction  as  has  been  shown  by  D.  Shah  et  al 
(2004). 

Most  of  the  RI  models  used  deal  with  monotonic  loading  response,  but  here  we  apply  our 
approaches  to  cyclic  loading  schemes  and  verify  our  simulation  predictions  against  test 
results  obtained  on  PWA  1484  Ni-based  superalloy  at  isothermal  conditions.  There  are 
several  well  known  observations  that  should  be  modeled  while  analyzing  the  cyclic  tests: 

(i)  Bauschinger  effect  -initial  reduction  of  the  reverse  proportional  limit  and 

(ii)  Cyclic  softening-  reduction  of  initially  high  hardening  rate  due  to  cyclic 
transient  effects  (Watanabe  et  al  2002,  Xu  and  Jiang  2004). 

The  Bauschinger  effect  is  controlled  by  the  evolution  equation  for  the  back  stress  co{t) 
(see  the  eq.  7).  During  the  reverse  cycle  the  back  stress  is  still  not  yet  steady  and  the 
absolute  value  of  the  effective  stress  r  -  co\  reaches  the  critical  value  earlier  than  during 

the  initial  “forward”  loading  step.  Experimental  results  (C.C  Lin  1995)  indicate  that  the 
cyclic  softening  effects  are  active  only  during  transient  regime  which  is  approximately 
equal  to  the  time  needed  for  the  back  stress  to  change  sign.  We  simulated  this  effect  by 
adding  a  latent  cyclic  softening  term  proportional  to  the  tenns  depended  on  sign{co  ■  ya  ) 
as  follows: 


Sa  =  j(2  -  r/ +  rjsign(coa  -ya)) 


(30) 


where  the  parameter  r)  has  to  be  adjusted  to  reflect  the  intensity  of  cyclic  softening.  In 
our  calculations,  we  used  rj  =  1 ,  which  means  that  there  is  a  non-hardening  plateau  during 
transient  unloading.  The  latent  hardening  law  has  been  modified  as  follows: 


sa  =  ga  h 


f  „<*\p 

1-^r 


V 


(31) 


J 


If  the  back  stress  and  the  effective  RSS  have  the  same  signs  the  dislocation 
micro  structure  deforms  in  a  similar  to  monotonic  loading  and  the  parameter  g  =  1 . 
Otherwise,  if  the  signs  of  the  back  stress  and  slip  shear  rate  are  opposite  the  dislocation 
cells  structure  is  dissolving  and  the  hardening  rate  slows  down.  It  is  important  to  note  that 
during  any  monotonic  loading  parameter^  =  1  and  Eq.  31  becomes  the  classical  Asaro 
hardening  rule  again. 


Results  of  the  model  implementation  against  test  data  are  shown  in  Figs  13-15.  Figure  13 
(a,  b,  c)  shows  a  comparison  of  model  prediction  against  experimental  data  for  <00 1> 
single  crystal  PWA1484  cyclically  loaded  up  to  1%  of  strain  amplitude  at  a  temperature 
of  870  C.  Due  to  strain  hardening  the  overall  loops  shift  up  in  the  stress-strain  plane  and 
the  hysteretic  loop  is  getting  narrower.  We  present  the  modeling  results  together  with 
measured  test  data  for  the  first  cycle  (Fig  13  a)  and  cycle  #10  (Fig  13  b),  as  well  as  cycle 
#  30  (Fig.  13  c)  where  the  cyclic  loop  is  stabilized.  One  can  see  that  simulation  results  are 
extremely  close  to  the  test  observations.  Fig.  14.  shows  stress-strain  relations  for  <1 1 1> 
oriented  single  crystal  cyclically  deformed  using  strain  control  to  the  mechanical  strain 
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range  of  0.8%  (actual  magnitude  is  0.41%)  at  temperature  of  870  C.  Due  to  the  high 
Young’s  modulus  the  yield  is  more  pronounced  along  the  <1 1 1>  orientation  and  the  loop 
is  wider  open  than  at  the  similar  loading  conditions  along  the  <00 1>  crystallographic 
direction.  Comparison  of  predicted  vs.  measured  data  demonstrates  high  quality  of  the 
model.  Modeling  results  obtained  for  <  1 23  >  single  slip  orientation  compared  against 
test  data  are  shown  in  Fig.  15  for  the  isothennal  conditions  of  870  C  and  strain  range  of 
0.8%.with  R=-l.  The  deformation  starts  with  (1 1 1)[01  1]  single  slip  system.  The  loop  is 
extremely  wide  due  to  relatively  low  yield  value  for  the  single  crystal  along  single  slip 
crystallographic  direction.  Error  of  the  numerical  prediction  is  less  than  10%  which  is 
good  for  a  single  slip  crystallographic  orientation.  Model  predictions  can  be  improved  by 
refining  tension-compression  asymmetry  properties  in  the  model.  We  did  not  pay 
attention  to  this  issue  in  this  work.  The  numerical  results  demonstrate  that  our  developed 
elastic-plastic  model  is  capable  of  obtaining  extremely  good  predictions  for  both  the 
monotonic  and  the  cyclic  plastic  response  of  a  LI2  single  crystal. 


6.  Concluding  Remarks 

The  constitutive  model  developed  has  been  implemented  in  the  commercial  finite 
element  software  ANSYS  as  a  material  user  routine  to  predict  yield  anisotropy  and  yield 
-thennal  dependence.  The  equations  governing  the  mechanical  response  have  been 
calibrated  using  existing  experimental  data.  The  model  predicts  the  crystallographic 
lattice  rotation  during  deformation,  which  is  important  during  material  processing  and 
cyclic  ratcheting,  especially  around  geometrical  features  such  as  cooling  holes  in  single 
crystal  turbine  blades.  Our  rate-independent  cyclic  crystal  plasticity  fonnulations  are 
designed  for  cyclic  and  non-isothermal  loading  conditions.  They  can  uniquely  determine 
the  amount  of  shear  along  active  slip  systems  at  each  increment.  The  approaches 
developed  are  numerically  robust  and  efficient,  allowing  numerical  analysis  of  tens  and 
even  hundreds  cycles  providing  a  working  tool  for  low  cycle  fatigue  (LCF)  and  thermal 
mechanical  fatigue  (TMF)  prediction  analysis. 

We  have  developed  a  rate  independent  model  that  can  be  readily  applied  to  TMF  over  an 
extremely  wide  range  of  conditions  which  naturally  reduces  to  our  SVD  approach  as  a 
limiting  condition. 

Based  on  our  modeling  results  and  their  comparison  with  experiments  it  is  possible  to 
conclude  the  following: 

(1)  The  major  deformation  mechanisms  of  high  temperature  creep  of  Ni-based  single 
crystal  superalloy  are  octahedral  (111}(110)  and  cube  {00 1 }( 1 10) 
crystallographic  slip. 

(2)  Single  slip  crystallographic  orientations  defonn  mostly  non-homogeneously, 
which  leads  to  necking  and  subsequently  shorter  life. 

(3)  Ovalization  of  the  plastically  deformed  single  crystal  specimen  might  be 
correlated  with  time  to  failure 

(4)  Strain  hardening  and  cyclic  softening  can  be  accurately  predicted. 
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The  combined  theoretical-numerical-experimental  study  of  single  crystal  PWA  1484  Ni- 
based  superalloy  reported  here  represents  the  steps  at  understanding  the  difficult  and 
important  topic  of  visco-plastic  deformation  in  LI2  systems;  it  holds  substantial  promise 
through  future  work  and  through  further  refinement.  In  particular,  the  model  will  be  used 
to  more  accurately  account  for  the  thermal-cyclic  effects  and  for  creep-plasticity 
interactions. 
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Table  1. 


Slip  systems  operative  in  PWA1484  at  high  temperature 

Octahedral 


a 

kJk*] 

1 

(111)[01 1] 

2 

(i  i  i)[Toi] 

3 

(iii)[0iT] 

4 

(i  I  T)[oT  l] 

5 

(i  T  T)[ToT] 

6 

(i  T  T)[i  10] 

7 

(T  l  T)[0i  l] 

8 

(TiT)[ioT] 

9 

(T  l  T)[T  1 0] 

10 

(TTi)[0TT] 

11 

(T I  i)[ioi] 

12 

(TTi)[T  10] 

Cube 


a 

13 

(100)[011] 

14 

(100)[01 1] 

15 

(010)[101] 

16 

(oiO)[ioT] 

17 

(001)[1 10] 

18 

(00i)[i  To] 

16 


Time  Independent  Plasticity  -  <001> 


Strain 


Figure  1 .  Typical  stress-strain  curves  obtained  with  rate-independent  model  with 
different  values  of  parameter  ko. 


Figure  2.  Variation  of  Elastic  Properties  with  Temperature. 


Figure  3.  Variation  of  Initial  slip  resistances  for  octahedral  and  for  cube  slip  systems  with 
Temperature. 


Figure  4.  Elements  of  disturbed  lattice  without  constraint  and  with  stiff  constraint  (from 
Dieter  1988)  and  two  types  of  boundary  conditions  applied  to  single  elements  to  simulate 
simple  tension. 
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Figure  5.  Stress-strain  relations  for  simple  tension  along  (a)  <00 1>  and  for  (b)  <1 1 1> 
crystallographic  orientation 


true  tension  strain 


Figure  6.  Stress-Strain  relations  for  simple  tension  along  <  236  >  crystallographic 
orientation  obtained  by  both  reported  methods  of  rate-independent  plasticity. 
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Figure  7.  Slip  systems  activity  for  simple  tension  of  single  crystal  PWA  1484 
along  <  236  >  crystallographic  orientation. 


Strain 

Figure  8.  Stress-Strain  relations  for  simple  tension  along  <  236  >  crystallographic 
numerically  orientation  obtained  by  applying  different  boundary  conditions  (single 
element  type  A  and  type  B)  and  for  an  element  in  the  3744  -element  cylinder  model. 
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Figure  9.  Crystal  lattice  re-orientation  of  the  single  LU  crystal  due  to 
simple  tension  along  <  236  >  crystallographic  axis;  (a)  Type  A 
boundary  conditions,  (b)  Type  B  loading  conditions,  and  (c)  central 
element  of  the  cylinder  model. 
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Figure  10.  Cylinder  FEM  used  for  the  simulation  of  single  crystal  tensile  sample 
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Figure  11.  Stereographic  view  of  the  necked  <  236  >  single  crystal 
specimen  after  overall  tension  to  20%. 


Figure  12.  Elliptical  cross-section  in  the  necked  part  of  the  specimen 


Stress  (MPa) 


4000- 


25 
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Figure  13.  Predicted  stress-strain  relations  against  experimental  data  for  strain-controlled 
cyclic  test  up  to  1%  strain  along  <00 1>  crystallographic  direction  at  870  C.  Results  are 
shown  for  the  first  cycle  (a),  cycle  10  (b),  and  cycle  30  (c). 


26 


Stress  (Mpa) 


Figure  14.  Predicted  stress-strain  relations  against  experimental  data  for  strain-controlled 
cyclic  test  up  to  0.8%  strain  range  along  <1 1 1>  crystallographic  direction  at  870  C. 
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Figure  15.  Predicted  stress-strain  relations  against  experimental  data  for  strain-controlled 
cyclic  test  up  to  0.8%  strain  range  along  <  1 23  >  crystallographic  direction  at  870  C. 
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